Mon. Not. R. Astron. Soc. 000,[T]{T6](2002) Printed 19 August 2010 (MN EGgX style file v2.2) 



Monte Carlo Simulations of Star Clusters - VII. The globular cluster 
47Tuc 



o 

W). 

< 
oo 



o 

6 

CO 

o3 



> 
OO 

O 

cn 

od 
O 

o 



Mirek Giersz 1 * and Douglas C. Heggie 2 

1 Nicolaus Copernicus Astronomical Centre, Polish Academy of Sciences, id. Bartycka 18, 00-716 Warsaw, Poland 

2 School of Mathematics and Maxwell Institute for Mathematical Sciences, University of Edinburgh, King's Buildings, Edinburgh EH9 3JZ, UK 



Accepted Received . . . ; in original form . 



ABSTRACT 

We describe Monte Carlo models for the dynamical evolution of the massive globular cluster 
47 Tuc (NGC 104). The code includes treatments of two-body relaxation, most kinds of three- 
and four-body interactions involving primordial binaries and those formed dynamically, the 
Galactic tide, and the internal evolution of both single and binary stars. We arrive at a set of 
initial parameters for the cluster which, after 12Gyr of evolution, gives a model with a fairly 
satisfactory match to surface brightness and density profiles, the velocity dispersion profile, 
the luminosity function in two fields, and the acceleration of pulsars. Our models appear to 
require a relatively steep initial mass function for stars above about turnoff, with an index 
of about 2.8 (where the Salpeter mass function has an index of 2.35), and a relatively fiat 
initial mass function (index about 0.4) for the lower main sequence. According to the model, 
the current mass is estimated at 0.9 million solar masses, of which about 34% consists of 
remnants. We find that primordial binaries are gradually taking over from mass loss by stellar 
evolution as the main dynamical driver of the core. Despite the high concentration of the 
cluster, core collapse will take at least another 20Gyr. 
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1 INTRODUCTION 

This is one of a series of papers in which we attempt to construct 
dynamical evolutionary models of old star clusters in the Milky 
Way. The fi rst in the series concerne d the massive, southern cluster 
lu Centauri dGiersz &~ Heggie 2003). That study was a pilot project 
which used a more primitive version of the code than is avail- 
able now, especially with regard to stellar evolution. Nevertheless it 
yielded quantitative theoretical results on mass segregation, which 
is usually neglected in studies of this cluster because of its very 
long relaxation time. Our updated code, more or less in its present 
form, was then tested in applica tion to the old open cluster M67 
dGiersz. Heggie, & Hurlevll2008l) . by comparison with existing N- 
body modelling. Then we returned to the globular clusters, for 
which iV-body modelling is eith er impossible or very lim ited. Our 
model of the nearby cluster M4 (He ggie & Gierszll2008al) showed 
that this object appears to have passed through core collapse, de- 
spite its uncollapsed (King-like) surface br ightness profile. Turnin g 
next to the rather similar object NGC6397 dGiersz & Heggij2009h . 
which does have a profile of the type thought to be typical of a 
post-collapse cluster, we found that it should be in a similar evo- 
lutionary phase to M4, and concluded that the difference between 
the two clusters was best explained as a result of fluctuations. De- 



E-mail: mig@camk.edu.pl (MG); d.c.heggie@ed.ac.uk (DCH) 



tailed TV-body modelling dHeggie & GierszH2009l) confirmed that 
a post-collapse cluster like NGC6397 should exhibit oscillations, 
on a time scale of order 10 8 yr, as a result of which it could some- 
times look like M4, and sometimes like NGC6397. In summary, 
each time we have looked at a specific globular cluster with our 
technique, we have discovered something slightly surprising: mass 
segregation in u Cen, the postcollapse status of M4, and core oscil- 
lations in NGC6397. 

In the present paper we turn our attention to the other famous, 
massive, southern globular cluster - 47 Tucanae (NGC104). Like 
M4 and NGC6397 there is a wealth of observations, even though 
it is considerably more distant. It is also much richer, reducing 
some statistical errors in such data as the surface brightness dis- 
tribution. It was selected for detailed study during discussion at 
the programme "Formation and Evolution of Globular Clusters", 
held at the Kavli Institute for Theoretical Physics, Santa Barbara, 
California, USA, in January 2009, just as M4 had been selected 
at the MODEST-5 meeting ("Modelling Dense Stellar Systems") 
at McMaster University, Hamilton, Canada, in August 2004. From 
the point of view of dynamical evolution, it differs from most of 
our earlier targets because of its long relaxation time (which is a 
consequence of its much greater mass). This makes it dynamically 
younger, raising the expectation that it should be easier to model. 
But the long relaxation time should also imply that its present state 
is more directly influenced by its initial conditions. 
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This paper first reviews both the observational data on 47 Tuc 
and earlier work on its dynamics (Sec. 2). Then in Sec. 3 we dis- 
cuss our techniques, in terms of the simulation of dynamical evo- 
lution, the calculation of mock observational data and our search 
for satisfactory initial conditions. In the following section, Sec.4, 
we present our optimal model, in the sense that, among the mod- 
els we have studied, it is representative of those that provide the 
most satisfactory fit to a wide selection of the available observa- 
tional data. Finally (Sec. 5) we consider such issues as the plausi- 
bility of our initial conditions, how these relate to current ideas on 
the earlier stages of cluster formation and evolution, and the mech- 
anisms which dominate the dynamical evolution of the cluster at 
the present day. Section 6 sums up. 



2 A REVIEW OF 47 TUC 



2.1 Dynamical models 



As with all globular star clusters, dynamica l models of 47 Tuc 
fall into two clas ses i Heggie & Giersz|[2008bl) : static models, such 
as King models dKindl 19661) , and dynamical evolutionary models, 
which are the focus of our own contribution. (We exclude here the 
practice of fitting a template to such data as the surface brightness 
profile. Though often referred to as King models, and useful as 
they are for the measure ment of such parameters as the core radius 
(e.g. lCalzettietalJ 19931) . they have no real dynamical basis.) In this 
subsection we survey existing models of these two kinds, though 
we often have to refer to various kinds of observational data on 
whi ch we expand in later subse c tions. 

Illlingworth & Illingworthl d 19761) fitted a single-mass King 
model to a composite surface brightness profile constructed from 
photometric data and star counts. Sim ple though i t is, s uch a 
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model was also used more recently by iFreire et al. feOOlh for a 
study of pulsar accelerations; and by iMcLaughlin et al.H2006h 
for the interpretation of the velocity dispersion profile, derived 
from proper motions and radial velocities. Several investigators 
have constr ucted multi-mass generalisa t ions of the King model 
for 47 Tuc foa Costa & Freernanl 1 19851; iMevlan & Mayo j[l9li ; 



iMapelli et al ll2004l) . lMevlarj dl988Ml989h and lReej d 1 996h added 
anisotropy, in the style of Michid dl963l) . 

A number of other, rather different static models have been 
constr ucted for 47 Tuc. It was pointed out by ISau Fa & Pedronl 
that a better fit to the surface brightness profile could be 
obtained by construction of a model with a slightly different choice 
of distribution function from that of King. The rotation of 47 Tuc 
has been the subject of considerable study, but almost all dynam- 
ical model s ignore i t. One exception is the pair of models con- 
structed bv lDavous Finally in this s ection we mention the 
non-parametric analysis of iGebhardt et al .1 dl995t) . which derives 
the distribution function from the observational data. 

More in line with the aims and methods of this pape r are 
the dynamical evo lutionary models of iMurphv et al.l dl998h and 
iBehler et al.l d2003l) . Though few details are available in these ab- 
stracts, the technique is a multi-mass Fokker-Planck code, with a 
discretised, piecewise power-law mass function and heavy rem- 
nants. The models are able to account for the stellar mass func- 
tion, star-counts, the velocity dispersion profile, pulsar accelera- 
tions, and the radial distribution of neutron stars. The models fit 
better than a King-Michie model. 



2.2 Observational Data 

Our purpose in this subsection is to review something of the wealth 
of observational data which exists for 47 Tuc. Not all of this data 
is compatible, and one of our aims is to explain where and why we 
have selected one dataset over another. 



2.2.7 Surface brightness profile and star counts 

In several of the foregoing investigations the spatial distribution 
of matter in 47 Tuc has been specified in terms of a surface 
brightness profile, and two major compilations exist dMevlanl 19881 ; 
ITrager etal.ll 19951). The later o f these is largely a superset of the 
former, though iMevlanl Jl988l) includes star count data in J from 
|Pa Costal d 19821) which are excluded fr om the later c ompil ation. 
The two compilations agree quite well; iTrager etal I d 19951) also 
provide an analytic fit to thei r data, and the rms difference between 
this and the data in iMevlar] dl988l) is 0.28 mag. A large fraction 
of this is contributed by star counts at large radius, including the 
aforementioned J data. 

Despite this apparently satisfactory state of affairs, some cau- 
tion must be exercised in the construction of a composite profile 
from such disparate data. Da Costa himself (his Fig. 7) shows that 
the slope of the density profile from different plates (differing es- 
sentially in the range of masses of the counted stars) is different, 
in a manner consistent with mass segregation. A further complica- 
tion i s the colour gradient exhibited by 47 Tuc dChun & Freemanl 
1 1979b . which has been traced to a radial variation in the fraction of 
light contributed by bright giants dFreemanlll985l) . One implication 
of these results is that it is not strictly possible to reconstruct the 
surface brightness from star counts by a constant shift, as is done in 
these compilations. Furthermore, beyond about 5' the only genuine 
p hotometric data in these co mpilations are the drift scan measures 
of iGascoigne & Burr! dl95^) , and an inspection of their Fig.l does 
not inspire much confidence in the results at these radii. With this 
consideration in mind, in our model ling we have not re lied solely 
on the surface brightness profile of ITrager et alj dl995h . but have 
also compared our models with counts f rom two sources, (i ) the 
two deepest V data in the c ompilation of Trager et al. I dl995h . i.e. 
plates 3407 and 3754 from |Pa Costal d 19821) . and 7ii) the surface 
density profile (based on star coun ts for stars above turnoff m ass in 
the innermost 82.5") published by IMcLaughlin et alj d2006l) . 

One significant de tail about our handli ng of the data in 
IMcLaughlin et al.ld200o1) concerns the passband. IMcLaughlin et alj 
1 20061) give surface densities for stars above about turnoff; more 
precisely, for st ars with mF475W < 17. 8, corresponding roughly 
to V = 17.65 dMcLaughlin et ak I l2006l , Sec.4. 1). These are em- 
pirical values, whereas our model creates V magnitudes based on 
model atmospheres. We have checked several synthetic databases, 
and find that mF47sw — V is considerably larger. For example for 
Padova isochrones (http://stev.oapd.inaf.it/cgi-bin/cmdl it is 0.284, 
for the bluest stars around tumoff at the appropriate age and metal- 
licity, and for the appropriate Dartmouth models it is still larger. In 
Sec. 14.11 we present results based on both the empirical and syn- 
thetic values. 



2.2.2 The velocity dispersion profile 

For velocities we initially relie d on the line-of-sight velocity disper- 
sion profile of IMevlanl dl988l) . which was based on a catalogue of 
169 radial velocities of giants. Nevertheless much larger data sets 
have been acquired since, and we begin this section by reviewing 
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Gebhardtetal(1995) 
Meylan(1988) 
Lane et al (2010) 




Figure 1. Line-of-si ght velocity dispersio n profiles of iMevlanl Jl988h . 
IGebhardt et"al] Il995h and lLane et al] 42O10t) . For the first of these the verti- 
cal error bars are as in the original publication, and the horiz ontal error bars 
give th e range of projected radii in which the stars fall. For Gebhardt et al. 
119951) we have copied curves in their figure which give estimates of the 
velocity disp ersion profile and its 90% confidence interval. For the most re- 
cent data set iLane et alj2010h we have copied the points and vertical error 
bars from their Fig 1 1 . 



more recent work, with a view to determining if the older data is 
adequate for our purpose. 

Considerable interest has cen tred on the subsequent discovery 
dMevlan. Dubath. & Mavorj|l99lh of two high- velocity stars with 
line-of-sight speeds (relative to the mean of the cluster) approx- 
imately 4 times larger than the central line-of-sight velocity dis- 
persion. It can be argued that there is nothing particula rly ab- 
normal about these dPetersoJl994l : lMcLaughlin et al.l2006l) . They 
do, however, have the effect of increa sing the central line-of-site 
veloc ity dispersion from 9.4±1.0km/s dMevlan. Dubath. & Mavorl 
Il99ll) to 11.6km/s. Therefore it would not be surprising if the cen- 
tral velocity dispersion of a successful model exceeds that of the 
older data. 

Radial velocitie s for 548 stars in the c entral region of 47 Tuc 
have been secured by[Gebhardt et "al] ( ll995l) . and the inferred veloc- 
ity dispersion profile is plotted along with the older data in Fig[TJ 
This figure confirms the result that the central vel ocity dispersion 
may be consi derably higher tha n in lMevla n Ul988h . 

Recentlv lLane et al] d2010h obtained a velocity dispersion pro- 
file based on t he largest spectroscopic sample ever o bserved for 47 
Tuc (see also lKiss et al]|2007l ; ISzekelv et alj|2007h . and we have 
included the results in the figure. Their results show that for the 
outer parts of the cluster the velocity dispersion appears to be ris- 
ing with radius, while the central value (obtained by fitting a Plum- 
mer model) is only about 9.6 km/s. Note, however, that one of 
the criteria for membership of 47 Tuc in their study was that the 
radial velocity must lie within a certain range of width about 45 
km/s (their Fig.l). T his would have excluded the two h igh- velocity 
stars discovered by dMeylan, Dubath. & Mavorl 1 1 99 ll) . and proba- 
bly suppresses the velocity dispersion at small radii. Twenty-five 
stars w hich, according to the criteria of membership in lLane et al.l 
are members of 47 Tuc, actually lie outside the conven- 
tional tidal radius. Several of these extratidal stars lie in sparsely 
populated parts of the colour-magnitude diagram (their Fig. 3), 
which suggests that they may not be members at all. It is possible 
that the membership criteria are insufficiently strict at large radii. 



On account of such arguments, for the work reported here we 
have relied mainly on the data of Gebhardt et al, and so we now dis- 
cuss it a little further. The sha rp drop in the Gebhard t et al data just 
inside 1 arcmin is interesting. IGebhardt et alj J 19951) regarded it as 
a sharp rise (as the radius decreases), and suggested tentatively that 
it may indicate a central pop ulation of dark objects. The analysis 
of lGebhardt & Fischerl i ll 9951) , however, implies that the M/L ratio 
inside this feature is little higher than at 5 arcmin. In fact it implied 
that there was a zone of very Jow M/L ratio around 1 arcmin (pre- 
sumably because the drop in velocity dispersion in this vicinity is 
almost Keplerian). Furthermore the 90% confidence limits on this 
result require an increase in M/L with increasing radius between 
about 1 and 4 arcmin. The lesson we draw from this discussion is 
that the features in the velocity dispersion profile may have little 
significance, and indeed the confidence limits themselves are sub- 
ject to some doubt. 

In our selection of models (Sec [3]l we have not used the wealth 
of data on proper moti ons dReesI 1 19961 : iKing & Anderson! 1200 ll ; 
iMcLaughlin et al.l2006l) . the bulk of which is confined to the inner- 
most 100". Nevertheless we shall compare some aspects of these 
data, such as anisotropy (Sec |5.2t . with one of our best models. 

The rotation of the cluster is evident in both proper motions 
and radial velocities, and has been studied in many of the above 
papers. Unfortunately, however, our Monte Carlo technique is re- 
stricted to the case of spherical symmetry and no rotation. We re- 
turn to this issue in Secl4.2l 



2.2.3 Luminosity and mass functions 

Ground-based data on the mass function dHesser et alii 19871) pro- 
vided an impo rtant constraint for the fitting of King-like models 
dMevlanl 19891) , and we shall see that luminosity functions continue 
to serve such a role in our work. Indeed the use of the luminos- 
ity function is preferable, as it places the comparison closer to the 
observational domain, and simplifies problems connected with the 
mass-luminosity relationship. We have, however, mainly used more 
recent d ata: the luminosity funct ions derived from HST observa- 
tions by [Anderson &King| l fl996T) at 1 and 14 core radii, where the 
core radius was taken to be 23" . The outer field goes down to mag- 
nitudes corresponding to stellar masses of about O.lAf©. Other in- 
formatio n on the luminosity function , though not used here, can be 
found in lde Marchi & Parescd dl9 95): Santi ag o . Elson. & Gilmorel 
d 19961) and lHowell. Guhathakurta. & Gillilandl fcOOCl). 



Hesseretal.ldl987h for 



We have retained the older data of 1 
discussion in Sec. 14.31 where we consider their "composite" lu- 
minosity function. As the authors themselves state, they did not 
attempt to correct for mass segregation when combining data se- 
cured at different radii, and it should not necessarily be regarded 
as a global luminosity function. Our interest in it is focused on the 
stars above turnoff, however. 



2.2.4 Binary fraction 

The binary fraction is very relevant for a variety of exotic objects, 
including blue stragglers and millisecond pulsars. Nevertheless it 
seems to be observationally rather poorly constrained in most glob- 
ular clusters, and 47 Tuc is no exception. The issue is complicated 
by mass segregation and the observational selection effects of the 
various techniques for detecting binaries. 

From inspect ions of col our-magnitude diagrams, 

Ide Marchi & Parescd d 1995b and Santiago, Elson. & Gilmord 
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d 19961) estimated that the fraction of "equal-mass" binaries in fields 
4.6' and 5' from the cluster centre was at least 5 %. Using the same 
basic approach, at the same distance of 4.6 , I Anderson! Jl997h 
found a fraction of at most 2% with mass ratio q > 0.5, in the 
magnitude range from V = 20.5 to 23.5. Again using broadly the 
same technique, Milone (pers.comm.) finds the fraction of binaries 
with q > 0.5 to be 0.015 ± 0.008, and the fraction of all binaries 
to be 0.026 ± 0.015. His measures refer to a region between 0.95 
and 2.4 arcmin from the centre, and he points out that they take no 
account of the recently discovere d intrinsic broadening of the main 
sequence jAndersonetalJl2009l) . In the central regions (< 90") 
lAlbrow et al. I J200lh estimated a binary fraction of about 14%, 
based on observations of eclipsing binaries and some modelling. 
The stars observed range in V fro m about 16 to about 23 , and 
it is known from modelling (e.g. iHeggie & Gierszl l2008ah that 
mass segrega t ion is particularly pronounced among bright objects. 
Knigge et alj u2008r> estimated a similar binary fraction among 
white dwarfs in the core. 

Even though these estimates refer to somewhat different pop- 
ulations of binaries, and in some cases at different locations, the 
binary fraction seems to be the least well constrained of the obser- 
vational quantities we have discussed. 



2.2.5 Pulsars 



3 THE SIMULATIONS 



3.1 Model description 



It is estimated dHeinke et ail 120051) that 47 Tuc contains about 
300 neutron stars. Am ong these are 23 observed millisecond pul- 
sars dFreire et al.ll2003lfl . Since our modelling focuses on dynam- 
ical issues, the interest of these objects is two-fold. First, some 
previous models of 47 Tuc have invoked such heavy remnants 
as b eing necessary to accou nt for the central velocity dispersion 
(e.g. lMevlan & Mavodl 19861) . Se cond, well-obs erved pulsars act as 
probes of the gravitational field dPhinne\ill992l) , in a way that we 
now summarise. 

Sixteen of the millisecond pulsars have timing solutions which 
include measurements of the logarithmic period derivative P/P, 
where P is the pulse period. One contribution to this is intrinsic, 
caused by the spin-down of the pulsar, but the size of this contri- 
bution (which is positive) is estima ted to be of order (P/P)i„t ~ 
0.5 x 10~ 17 /sec dFreire et alj|200ll) . In addition there is a Doppler 
contribution from the relative acceleration of the pulsar and the ob- 
server, and the main contribution is generally that due to the matter 
in the cluster itself. For a spherically symmetric cluster the magni- 
tude of this term is GM(r) cos 8/{cr 2 ), where M(r) is the cluster 
mass within the radius r of the pulsar, G, c are the usual physi- 
cal constants, and 9 is the angle between the line of sight to the 
pulsar and the outward radial direction at its location in the clus- 
ter. Thus the spin period of a pulsar located on the far side of the 
cluster appears to be decreasing, unless the contribution from in- 
trinsic spin-down is too large. Of course only the projected radius 
r p = r sin 8 is known with any accuracy. For a given cluster model 
and given r p , however, the extreme values of the Doppler contribu- 
tion can be calculated. Then, all observed values of P/P should lie 
between these extrema, except for the small intrinsic contribution. 
This is illustrated for one of our models in Figll II 
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As in previous papers in this series, the code we use i s an updated 
versio n of the Monte Carlo code developed in lGierszl J 1 998L l200ll . 
1 20061) . It includes synthetic stellar e volution of single and binary 
star s using prescriptions describ ed bv lHurlev. Pols. & Toutl l l2000l) 
and iHurlev. Tout. & Pols! d2002h - The only relevant updates to re- 
port here are the inclusion of (i) exchange interactions (through the 
use of cross sections), and (ii) natal kicks for both neutron stars 
and black holes (which had been confined to single neutron stars in 
previous versions of the code). 

For initial conditions we adopt the assumptions given in Table 
[T] The table also includes an indication of the range of numerical 
values which were explored in the search for a model which is de- 
scribed in Sec l3.3l Some of our choices require a little elaboration 
here. 

(i) The initial mass function (a two-part power law) depends 
on five parameters, but only four were adjusted in the search for 
a model. For the low-resolution initial search (Sec l3,3b the mini- 
mum mass was held at mi = O.IMq, as in previous papers in this 
series. It was kept at mi = O.OSA/q for the full-resolution models 
discussed in Sec|4] 

(ii) As in our work on M4 dHeggie & Gier sz 2008a), the binary 
parameters follow the prescription of iKroupal 1 1995b . except that 
the initial component masses are derived from the same initial mass 
function as the single stars in our model; its parameters normally 
differ in value from those in the mass function used by Kroupa. 
The component masses, period and eccentricity are then subject to 
"eigenevolution" and "feeding" algorithms. These have a physical 
basis and are designed to build in a number of correlations consis- 
tent with observational data. 

(iii) For the initial spatial distribution we have adopted King 
models, but models which (in general) underfill the initial tidal ra- 
dius. Thus there is a distinction between the tidal radius (deter- 
mined by the strength of the tidal field and the mass of the cluster) 
and the edge radius of the King model. We have also experimented 
with several other types of model, as described briefly in Sec |4.2| 

(iv) A constant tidal field strength is used, as i f the clus- 
ter were in a circular orbit. F or 47 Tuc iTucholke] d 19921) and 
iDinescu. Girard. & van Altenal d 19991) estimated an eccentricity 
about 0.17, i.e. the Galactocentric distance varies by approximately 
this fractional amount on either side of its mean value. Therefore 
our assumption has some approximate validity. 

(v) Though it is not strictly an issue of modelling, for the pur- 
poses of comparison between model and observations we need to 
choose values for the age and distance of the cluster. The adopted 
distance is the value given in Table [2] but for the main model de- 
scribed in Sec|4]we found it slightly advantageous to choose an age 
at the upper end of the range in the Table, i.e. 12 Gyr (see Sec|4~] 



3.2 Generation of observational data 



We have already discussed (Sec l2.2t the main kinds of observa- 
tional data which we shall attempt to fit with our model. Here we 
give the recipes used for the generation of the corresponding data 
from our model (and some additional data such as the core radius). 

At any time, the output of the Monte Carlo model consists of 
a list of data for each star or binary. This data includes the radius r 
(i.e. distance from the cluster centre), radial and tangential velocity 
(v r and vt), mass, ^-luminosity Ly in solar units, and colour of the 
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Table 1. Choice of initial conditions etc 



Property 


Distribution adopted 


Parameter 


Meaning 


Range (small- 
scale models) 


Range (full- 
scale models) 


Notes 


Number of stars 




N 




a 


1.8 - 2.2 X 10 K f> 




Structure 


Kins (1966) 


Wo 


Central potential 


4-11 


7-10 


c 


Tide 


Giersz et al (2008) 


n 


Tidal radius 


a 


72 - 116 pc 


See text Sec|3.Uiv) and d 


Radius parameter 




n/r h 




25 - 500 


21 - 606 


Determines (given rt) 


Mass function 


Two-part power law 


ni> 


Maximum mass 


2 - 150M Q 


2 - 150M Q 


KrouDa (2008) but with 


(single stars) 




nib 


Break mass 


0.86 - 0.9M Q 


0.75 - 1M Q 


different parameter values 






mi 


Minimum mass 


0.1M© 


O.O8M 








a 2 


Upper index 


2.3 - 5.3 


2.3-4.5 


Salpeter is 2.35 






ai 


Lower index 


0.4 - 1.3 


0.4 - 1.3 




Binaries 


Kiouna (1995) 


h 


Binary fraction 


- 0.20 


- 0.06 


See text Secl3.1lii) 


Mass segregation 


None 










But see the end of Sec 14. 21 


Natal kicks 


Gaussian 


a 


1 -dimensional dispersion 


190 km/s 


130 - 190 km/s 


neutron stars and black holes 


Age 








11 Gyr 


11 - 12 Gyr 




Metallicity 




Z 




0.0035 


0.003 - 0.0035 





Notes: a See text Sec |3.3l b Larger values were tried for truncated polytropic models, c Also lWoollevI (T954) and polytropic models (Sec |4.2> . d A wider range 
was tried for truncated polytropic models. 



star. For binaries the same data is available for both components, 
along with the semi-major axis. 

3.2.1 Surface brightness 

To construct the surface brightness at a given projected radius d, 
we replace each star by a uniform transparent spherical shell of 
radius r. Its surface (projected) luminosity density is given by 

Ey = — ^— T if r > d. The total surface luminosity 

27rH y r 2 - d 2 

density is the sum over all shells of radius r > d (which we also 
denote by Ey in what follows), and then the surface brightness is 
given (in V magnitudes per square arcsec) by 



M v = V Q -2.51og 10 E 



6480002 



(1) 



where V© is the absolute V-magnitude of the sun and Av is the 
visual extinction to the cluster. 

This method has a disadvantage that the brightness is infinite 
along a line of sight such that d = r for some shell. It would be 
possible to construct a realisation of the model by selecting the full 
position of each star at random on the corresponding shell, and then 
determine the surface brightness in a manner which is more akin 
to the observational approach. We have verified that this method 
produces similar results to the method described in the previous 
paragraph, and prefer that method (i.e. eq.lQJ), as it has the merit of 
introducing no additional source of noise. 

For the observational core radius we determine the smallest 
value of d such that the surface brightness is half that at the centre. 
(We also refer to a dynamical core radius in Sec |5.5l ) 



3.2.2 Surface density profile, binary fraction and luminosity 
function 

Our calculation of the surface density (i.e. number of stars per unit 
area on the sky) uses the same formulae as for surface brightness, 
as just described, except that Lv is omitted, and the conversion is 
to number per square arcsec and not V magnitudes. For the (local) 
binary fraction we compute the ratio nt/(nt + n s ), where nt is 
the projected density of binaries, and n s is the projected density 
of single stars. (Note, however, that the data in Table [4] give the 



global binary fraction.) The luminosity function is calculated in the 
same way as the surface density, but separates stars according to 
the appropriate magnitude bins. 

3.2.3 Line-of-sight profiles of velocity dispersion and proper 
motion 

Consider the point on a shell of radius r at which it is pierced by the 
line of sight of projected radius d < r. We assume the transverse 
velocity of the star is randomly orientated in the tangent plane to the 
shell, and then it is easy to show that the mean square line-of-sight 



velocity is given by v r 



1 2d 2 



+ — v t — 5-. For each shell within 



2 " 

the line of sight, this expression is then weighted by the surface 
density of the star, i.e. by a factor rj \/r 2 — d 2 , and the weighted 
average gives the line-of-sight velocity dispersion. 

For the line-of-sight velocity dispersion, stars fainter than ab- 
solute magnitude V — 6 are excluded from the sum. For proper 
motions similar formulae are used, with the appropriate magnitude 
limit (FigEJl. 



3.2.4 Pulsar acceleration 

For each line of sight, the line-of-sight component of the gravita- 
tional acceleration is computed at the location of each shell, and the 
maximum value found. 



3.3 Finding a model 

The starting point of our search for a model was a consideration of 
present estimates of the mass and dimensions of the cluster (Table 
O. The half-mass relaxation time is long, which suggests that the 
cluster has not lost a v ery large fraction of its initial mass. Indeed 
the formulae given by iBaumgard t & Making d2003h suggest that, 
even with its present mass, the dissolution time is of order 10 11 yr. 
This in turn would suggest that most of the mass lost so far by 47 
Tuc results from stellar evolution, and that its initial mass was less 
than twice its present mass. For much the same reason we began 
by assuming that its initial half-mass radius was not very different 
from its present value. Only in the core is the relaxation time short 
enough to suggest extensive evolution by relaxation. 
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Table 2. Data on 47 Tuc 



Rac 1 


7.4 kpc 


He 1 


4.5 kpc 


Eb-v 1 


0.04 


{Fe/Hf 


-0.76 


Total mass A/ 2 


1.1 x 10 6 M Q 


Half-mass relaxation time t r h 


3.0 x 10 9 yr 


Central relaxation time t rc 1 


0.9 x 10 8 yr 


Age 3 


x 10 9 yr 


,1 

' C 


0.40' 


1 


2.79' 




42.86' 



Sources and notes: ^Harris! Il99(j) . updated on 

i http://www.ph ysi cs.mcmaster.ca/~har ris/mwgc.dat on 19 January 2010; 2 
Meylan ( 1989): :i Gratton et al. (2003). though their results have been sim- 
plified slightly 

Guided by these ideas, we began the search for suitable initial 
conditions using small-scale models with N* = 10 4 or 4 x 10 4 
stars. These were scaled to have the sam e relaxation time as a full - 
scale model with N stars, as described in iHeggie & Gierszlfe008ah . 
Briefly, this is arranged by relating the radial scales of the two mod- 
els by 

r* _ ( n_ y/V iog(7Ag V /3 

R \Nj ^log( 7 JV) J ' ■ ' 

where 7 = 0.02. dPortegies Zwart et~al ] Jl999h used a similar trick, 
but with 7=1.) For iV* = 10 4 and N = 2x 10 6 (which is typical 
of the models that we eventually prefer) the ratio is R+/R ~ 3.68. 

In these small-scale models the age was taken as 1 1 Gyr, as 
in Table [2] and the minimum mass was taken to be O.IMq. These 
values differ from the values of 12 Gyr and 0.08A/© adopted for 
the main model in Sec[4] but the resulting differences are negligible 
for the purpose of the explorations discussed here. These models 
can be scaled to any value of N, which is not an input parameter 
of the model and is therefore not given in TableQ] Similarly we do 
not specify the range of r t for these models, as this also depends 
on the chosen scaling. 

Our first finding was that the canonical values for the initial 
mass function did not provide a satisfactory fit to the luminosity 
functions which we were trying to match (see Sec |2,2.3l >. We found 
that all our subsequent modelling of the luminosity functions (at 
least, below the turnoff) was satisfactory if we chose values closer 
to ai = 0.5 and = 0.9A/©, where these parameters are defined 
in TableQ] (Nevertheless, slightly different values were eventu- 
ally adopted for the model presented in Sec|4j see also Sec l4.5l ) 
Thus the low-mas s initial mass fu nction is flatter than in the canon- 
ical formulae of lKroup3 d2008l) . Strictly, only the value of Qi is 
approximately determined by this comparison; all we can say of 
mi, is that it is not much below the present-day turnoff mass, and 
the value chosen merely represents approximately the most modest 
change from the canonical values. 

Our conclusion about ai seems quite robust, and is illustrated 
in Fig(2] This shows a scaled version of a model close to the one 
to be presented in Sec[4] and one with identical parameters (see the 
caption to the figure and Table[3} except for the slope of the lower 
part of the IMF. Though a better fit could certainly be obtained 
by varying more than just one parameter, the figure does serve to 
illustrate the effect of a mismatch in the choice of IMF. 

In order to give a flavour of how our search for a model pro- 
ceeded from this point, we now explore a coarse grid of the main 
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Index 0.5: outer field + 
inner field x 

Index 1 .3: outer field * 
inner field □ 
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V magnitude 

Figur e 2. Observed luminosity functions at two radii (from 
lAnderson & Kindll996l) compared with two models. One is a scaled ver- 
sion of the model close to that presented in Sec|4] scaled to TV* = 40000, 
and the other is an identical model except for the power law inde x of the 
lower initial mass function, which is the "canonical" value of 1.3 iKroupal 
|2008|) . The actual initial parameters for the first model are given in Table[3] 

initial structural parameters, namely the initial half-m ass radius and 
concentration. (We adopt the models o f lKind Jl966h for the initial 
structure.) All the models in this particular discussion use a rela- 
tively large initial value of the tidal radius (r* = 500 pc, where 
we give the value for the small-scale model.) It should also be said 
that they represent a tiny subset of all the models we explored (see 
TableQ}. 

Some results of this representative mini-survey are sum- 
marised in Fig. [5] The results of each model were scaled (by the 
particle number) to fit the central surface brightness of 47 Tuc, and 
then this figure gives the line-of-sight velocity dispersion at a pro- 
jected radius of 1 arcmin. While these suggest that there is a range 
of relatively compact initial conditions (r ; * in the range from 1.5 
to 2 pc) yielding approximately the correct velocity dispersion (see 
FigQ}, further inspection reveals serious problems. 

Figs|4] and [5] compare the entire surface brightness and ve- 
locity dispersion profiles for the model with initial concentration 
Wo = 9 and scaled initial half-mass radius 2pc. It has been scaled 
to a system with N = 2 x 10 6 stars. While the central surface 
brightness and velocity dispersion at 1 arcmin match reasonably 
well, as required, the remainder of the profiles do not match 47 
Tuc at all. The surface brightness of the halo is too low, and th e 
same result is obtained if one checks star counts jPa Costa|[l982h . 
The model does, however, have nearly the correct core radius. The 
poor match of the velocity dispersion profile is even more reveal- 
ing. The deficiency at large radii may be attributable to the under- 
massive halo, as shown by the surface brightness profile. The rise 
in the velocity dispersion well inside the core radius, however, is 
in complete contrast with the observational result. This reveals the 
presence in the model of a segregated population of massive ob- 
jects. 

As an aside, we consider briefly the nature of this population. 
Because of the scaling, the model has a low escape velocity, though 
the velocity dispersion of natal kicks (of black holes and neutron 
stars) was unsealed, and so, not surprisingly, the population of neu- 
tron stars is very small. The compact central population consists 
almost entirely of white dwarfs, which make up 49% of the entire 
mass of the cluster. There is nothing abnormal about this, especially 
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Initial half-mass radius r h (0) (pc) 

Figure 3. Line-of-sight velocity dispersion at 1 arcmin for a number of 
scaled models. Different curves correspond to different initial concentra- 
tions. The abscissa is the initial half-mass radius of a small-scale model 
with N* = 10 4 stars, and is related to that of a full-scale model with N 
stars by eq.|2). The models have been scaled to the correct central surface 
brightness for 47 Tuc, and so the value of N differs from one model to 
another. 



_ 20 

E 



Scaled model 
Gebhardtetal(1995) 




Projected radius (arcmin) 



Figure 5. Line-of-sight velocity dispe rsion profile of a scale d model com- 
pared with the observational results of Gebhardt et al. 1 1995 ). The spike in 
the model is attributable to a single star at small radius. Only bright stars, 
with V < 5, were included. 




Scaled model 
Trageretal(1995) 



10 100 
Projected radius (arcsec) 



Figure 4. Surface b rightness profile of a scaled model compared with the 
composite profile of lTrager et alj Jl995h - 



given the relatively flat low-mass initial mass function, which sub- 
stantially decreases the mass of the main sequence below turnoff. 
In our favoured model, however (see Tabled, the white dwarf mass 
fraction is only 34%, and the central velocity dispersion agrees bet- 
ter (Fig[8}. 

Returning to the problem of finding a suitable scaled model, 
we can also view the results of this mini-survey from a different 
angle; and it turns out to be a more fruitful approach. Instead of 
scaling to match the central surface brightness, consider the re- 
sult if instead one scales to match the central line-of-sight veloc- 
ity dispersion. Then it is found that the resulting models are too 
faint centrally, and the (observational) core radius is too large, no 
matter how compact the initial conditions are. A better fit could 
therefore be obtained if expansion of the core could be suppressed. 
Various mechanisms are known to cause expansion of the core, in- 
cluding loss of mass by stellar evolution, binary heating, and the 



action of black holes. (The pa pers bv lChernoff & Weinberg! d l990h . 
iHeggie. Trenti. & Hutl | |2006|) . and iMerritt et alj | |2004|) . respec- 
tively, may serve as samples of the extensive literature on these pro- 
cesses.) We judge that the initial binary fraction adopted here (2%) 
is too low (for our choice of initial period and mass distributions) 
to cause significant core evolution; in any case the observational 
constraints (Sec |2.2.4"l > do not suggest that it can be much lower, 
and so we have not attempted to reduce the binary fraction. The 
remaining mechanisms for core evolution, however, all involve the 
more massive stars in the initial mass function, and we have con- 
sidered modifications here. 

Despite its implausibility, perhaps, our efforts initially focused 
on truncation of the high-mass IMF, by severely reducing mi ( Ta- 
ble QJ. Indeed values around 2M@ gave rather satisfactory scaled 
models, without any change in the slope of the upper mass function 
from the canonical value of 2.3. One problem with such a mass 
function is that no neutron stars are formed, and so we explored 
different ways of altering the IMF. Though any significant increase 
in the break mass always gave poorer results, results were almost 
as good if we raised the maximum mass again but compensated by 
steepening the power law index of the upper mass function. 

With only small changes in the remaining parameters, the pa- 
rameter values arrived at in this way led to the small-scale model 
specified in Table [3] Then, with further changes and much exper- 
imentation, we arrived at the full-scale model described in the 
next section. We summarise in Table Q] the entire ranges of all 
parameters which we explored, both with small-scale models (i.e. 
as described in this section) and full-scale models (as in Sec[4]l. 
We do not intend to suggest that all possible joint variations in 
these parameters were investigated; very often the extreme values 
in these ranges were tried only for one or very few choices of the 
remaining parameters. In principle it would be desirable to adopt a 
more systematic way of searching this ra ther large parameter space. 
iHarfst, Portegies Zwart. & Stoltd J^oTjgj) describe an iV-body study 
of the Arches cluster which illustrates such an approach, at least 
when adjusting a pair of parameters (in their case the initial virial 
radius and particle number.) 
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Table 3. An approximate small-scale model 



Parameter 



Number of stars 
Total mass (Mq) 
Tidal radius (pc) 
Half-mass radius (pc) 
Central potential of King model (Wo) 
Binary fraction 
Upper mass function index 
Lower mass function index 
Minimum, break and upper mass (Mq) 0.1, 0.9, 50.0 

Note: corresponding values for a full-scale model, where these are different, 
are given in brackets 



Initial value 

N* = 40000 (N = 2 x 10 s ) 
2.93 x 10 4 (1.47 x 10 6 ) 
295 (109) 
6.6 (2.4) 
8 

0.02 
3.3 
0.5 



Table 4. Initial and present-day conditions for model A 

t = 



IV 1 2.00 x 10" 

Mass(M ) 1.64 X 10 6 

Tidal radius (pc) 86 

Half-mass radius (pc) 1.91 
Half-mass relaxation time (Gyr) 0.70 

W 7.5 

Observational core radius (pc) 0.42 

Binary fraction 0.0220 

Index of upper IMF 2.8 

Index of lower (I)MF 0.40 

Minimum mass (Mq) 0.08 

Break mass (M Q ) 0.8 

Maximum mass (Mq ) 50.0 

White dwarfs (M ) 

Neutron stars (Mq ) 



t = 12Gyr 



1.85 x 10 H 
0.90 x 10 6 
70 
4.96 
3.66 

0.55 
0.0084 

0.36 2 



3.1 x 10 5 
313 



Notes: 1 number of single stars plus number of binaries; 2 measured best fit 
between 0.3 and 0.8Mq, main sequence stars only. 

4 A MODEL OF 47 TUC 

Unlike the models discussed in the previous section, the model we 
are about to present is a full-scale model of 47 Tuc. It will be seen 
that the profiles and other data are considerably smoother than for 
the models discussed hitherto. This standard model will be referred 
to as Model A. Details of its initial specification, along with a sum- 
mary of conditions at 12Gyr, are given in Table|4] Though the slope 
of the upper mass function may seem steep, st i ll high er values (in 
the range 3.75—4.5) were inferred by iMevlanl il989h . His result 
was based on the need to produce the appropriate mass of white 
dwarfs to account for the velocity dispersion profile. Our motiva- 
tion is rather different, as we have already mentioned at the end 
of the last section. Meylan, incidentally, took a much more nearly 
canonical value for the power law index of the lower mass function 
than ours (he took 1.2), but a similar break mass (he took 0.88Mq). 

It may be of interest to note that each full-sized model such 
as this takes less than a day on a single Dual-Core AMD Opteron 
2214 at 2.2GHz. 



4.1 Surface brightness and density profiles 

As can be seen in Fig(6j model A is somewhat bright at small 
radii. Quantitatively, ou r central value is ab out 14.1 in the units of 
the fig ure, while fits bv lTrager et al.l i 19951) and lMcLaughlin et al.l 
d2006l) give values of 14.42 and 14.26 respectively. Less evident, 
perhaps, because of the steepness of the profile at large radii, is 




Model A 
Trageretal(1995) 




10 100 1000 

Projected radius (arcsec) 



Figure 6. Surface brightness profil e of Model A compared with the com- 
posite profile of lTrageret afUl995t) . 



the fact that the model is too faint there. The mismatch approaches 
about 1 magnitude at the largest radii plotted, (though this is less 
severe than for the scaled survey model shown in Fig|4j and can 
be at least partially attributed to our treatment of the tide (see text 
below). 

Star counts offer an alternative approach to the comparison of 
the spatial structure of the model with 47 Tuc (FigO. The model 
is somewhat overdense at small radii. If the cutoff in V of the star 
counts is determined from synthetic model data, as described at 
the end of Sec l2.2.l1 the mismatch is approximately 0. 1 dex. This 
is approximately compatible with the excess brightness discussed 
above. 

At large density, however, the evidence from surface density 
gives a different impression from the surface brightness, as the fit 
with the observational data is now much better. Though the fit is 
still imperfect, it was evidence like this which for us confirmed the 
possibility that the composite surface brightness profiles might not 
be very reliable at large radii. We had been alerted to this by finding 
models which produced satisfactory mass functions (to judge by lu- 
minosity functions below the luminosity of the turn-off) but which 
seemed faint at large radii when judged by the surface brightness 
profile. 

Our treatment of the tide dGiersz. Heggie. & Hurley|2008l) is a 
further complication, as it implies that the most distant star actually 
lies some distance inside the tidal radius. For this reason, the fact 
that the surface density distribution of the model still falls below 
the observations at large radii is not necessarily a significant issue. 



4.2 Velocity dispersion profile 

This, presented in Fig(8] is perhaps the most unsatisfactory compar- 
ison for this model. Whatever the formal statistical result (by com- 
pariso n with the 90% confidence band estimated bv lGebhardt et al.l 
dl995h ). and however much it improves on the model shown in 
FiglS] it still gives the qualitative impression of not conforming to 
the shape implied by the observations. This is particularly notice- 
able in three places: at the smallest radii, where it might seem too 
high, and around 2' and 10', where the model velocity dispersion 
seems too high and too low, respectively. Besides possible short- 
comings of the model, a number of relevant factors should be borne 
in mind, and we consider them in turn. 
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Da Costa (1982) 
Model A (empirical calibration) 
Model A (synthetic calibration) 
McLaughlin et al (2006) 
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Figure 7. Surf ace densit y profi le of Model A compared with some of the 
star counts of |Pa Costal jl982l) . which have a limiting apparent magni- 
tude of my = 19.45 , and the surface density of stars above turnoff 
from McLaughlin et al. (2006). Details of the "calibration" are given in 
Secs l2.2.T1 and l4. ll For the model, we switch from one limiting magnitude 
to the other at 100". 



Looking first at small radii, up to about 0.5', we see that any 
well-fitting model with a relatively flat (or, at least, non-increasing) 
velocity dispersion profile must pass close to the 90% confidence 
limits at very small radii and/or around 0.4'. Therefore the fact that 
our Model A does so does not cast doubt on its correctness. 

The second point to note is that, as mentioned in Sec |2.2.2l 
our model takes no account of the rotation of the cluster. The 
maximum mean line-of-sight (rotational ) velocity, about 6 . 5km/s , 
occurs at a radius of about 5- 6 arcmin dMevlan & Mavod[l986h : 
the Fabry-Perot observations o f lGebhardtetal.Ul995l their Fig.l 1) 
extend out to slightly smaller radii, but they give a value above 
7.5km/s already at about 4 arcmin. Thus the rotation is large in the 
second of the three regions where the mismatch seems most pro- 
nounced. The influence of rotation may also extend to larger radii: 
it is worth no ting that the singl e mass, rotating Fokker-Planck clus- 
ter models o f lKimetalj|200a, their Figs 1,11 and 12) clearly show 
that rapidly rotating clusters evolve faster and have larger half-mass 
radii than non rotating clusters. This will put more mass into the 
outer parts of a cluster, and increase the velocity dispersion there. 

The third factor to bear in mind is the influence of the 
tide at large radii. We treat its action as a cutoff, which is 
an increasingly rough approximat ion as the tidal boun dary is 
approached. (The observations of iGebhardt et al. I i 19951) extend 
to almost half th e observed tidal rad ius.) It has been known 
since the work of iDrukier et al. I dl998h that the velocity disper- 
sion profile in some clusters flattens off towards the tidal ra- 
dius, presumably because of the effects of the galactic tide. Such 
an effect is cer t ainly observable in A^-body simulations (e.g. 
Giersz & Heggid 1 19971; ICapuzzo Dolcetta. Pi Matteo. & Miocchil 
20051 : iKupper et all 20 id) an d has also been confirmed recently for 



47 Tucbv lLaneetalJ l l2010l) (see Fig[TJ Therefore the apparently 
excessively rapid decline of the velocity dispersion at large radii 
in the model may be an artefact of our simple tidal treatment. We 
note in t his res pect that the multi-mass anisotropic King model of 
iMevlanl Jl988h did a much better job in fitting the outer velocity 
dispersion and its trend. The very large half-mass radius in that 
model (Sec|5~TTl is no doubt a factor, placing more mass at large 
radii and elevating the velocity dispersion there. The fact is, how- 




Projected radius (arcmin) 



Figure 8. Line-of-sight velocity dispersion profile of Model A c ompared 
with the 90% confidence band estimated bv lGebhardt et alj ll995t) . 



ever, that our models, which start as compact King models, do not 
evolve as far as those constructed by Meylan. Also, it seems likely 
that some of the velocity dispersion which Meylan's model fitted 
is due to processes which were not included in that model, i.e. the 
tidal field. 

It is worth reporting here that we have attempted, with very 
limited success, to improve the outer velocity dispersion by con- 
structing models which might enhance the mass density at large 
radii. In particular we tried 

(i) Woolley models dWoollevlll954l) . which employ a truncated 
isothermal distribution, in contrast to the lowered isothermal distri- 
bution of the King model; their concentration can be specified in 
much the same way, i.e. by the scaled central potential Wo; 

(ii) polytropic models, with polytropic index exceeding the 
Plummer value of 5; the models were truncated at the tidal radius 
and revirialised globally, just as would be done with a Plummer 
model; and 

(iii) initially mass segregated m odels generated according to ei - 
ther (a) the prescription given bv iMcMillan & Vesperinil d2007h. 
with time delay of stellar evolution up to 2 5 My r; or (b) the 
method of iBaumgardt. De Marchi. & Kroupal d2008h . with mass 
segregation down to 2 Mq . 

Mass segregated models do differ in the central part of the system, 
generally slightly lowering the velocity dispersion, but in the region 
farther out than about 5' they are indistinguishable from models 
without primordial mass segregation. Polytropic models in partic- 
ular provided a slightly better fit to the velocity dispersion profile, 
though the improvement is too marginal to justify abandoning King 
initial conditions. 



4.3 Luminosity and mass functions 

Luminosity functions in fiel ds at two radii are shown in Fig[9] in 
comparison with data from lAnderson & Kind dl99ot) . According 
to those authors, turnoff corresponds to about V — 4.1 on the scale 
of this diagram. Below turnoff, then, the agreement is fairly satis- 
factory in both fields, in terms of both overall shape and normalisa- 
tion. The quality of the normalisation is relevant to the discussion 
in Sec |4. 11 For instance the inner field is at a radius of about 23" 
where, as we have seen, the evidence of both number counts (for 
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Figure 9. Lumi nosity functions at two rad ii, in Model A and from the obser- 
vational data of lAnderson & Kinel ll996l) . Note that these authors specify 
the radii in terms of the core radius (1 and 14r c ), and we have assumed 
r c = 23". 



stars above turnoff) and surface brightness is that the model is too 
dense or bright by about 0.1 dex. Nevertheless the mean difference 
in the logarithm of the luminosity function in Fig(9]between V = 4 
and 8 is 0.02 ± 0.03 (standard error), i.e. the model is not signifi- 
cantly denser than the observations. 

This discussion is co mplicated by the fact that the data 
of lAnderson & Kind i ll 9961) do not extend much above turnoff, 
whereas the surface density data discussed in Sec l4.1l is confined 
to stars above turnoff, and they also dominate the surface bright- 
ness. To explore the luminosity function above turnoff, in FigllOl 
we compare our glo bal luminosity funct ion with the composite lu- 
minosity function of lHesser et al.l Jl987l) . We have somewhat arbi- 
trarily normalised the latter (by eye) in the upper main sequence, 
where we know from the evidence of Fig|9]that the model appears 
to be satisfactory. A satisfactory match continues up to about the 
prominent peak (at the location of the horizontal branch), but no 
further. Actually, the distribution along the giant branch is known 
to be a sensitive test of stellar evolution dBergbusch & VandenBerd 
1200 lh . compl icated in 47 Tu c by the spatia l variation of the giant 
branch itself dFreeman|[l98l ; iBailvnll 1994) . Though we therefore 
make no attempt to resolve the mismatches in this brightest part 
of the luminosity function, it plays a significant role in the surface 
brightness, and problems here may well contribute to the fact that 
the surface brightness and density of the model appear to be exces- 
sive. 

Let us return again to a consideration of stars below turnoff. 
The low-mass mass function is often discussed in terms of the 
power-law index over some convenient mass range. We have pre- 
ferred to carry out the comparison in the observational domain, i.e. 
in terms of the luminosity function, partly because the mass range 
over which the power law is fitted varies from one author to an- 
other, and partly because it is not always clear how the fit is ac- 
tually carried out. In any event, such comparisons depend on the 
choice of mass-luminosity relationship. Nevertheless we give one 
example, wh i ch is a comparison with the ground-based result of 
iHesser et alj Jl987h for the global mass function. They found an 
index of 1.2 ± 0.3 (Salpeter = 2.35) for masses corresponding to 
5 5} V 9, which corresponds to masses between about 0.5 and 
0.8Mq. We find at least that an index of 0.9 provides a reasonable 



Figu re 10. The global lum inosity function of model A in comparison with 
that o f lHesseret afll 19871) , The latter has been normalised to give a reason- 
able fit in the upper main sequence. 

fit in this range. Nevertheless, a smaller value provides a better fit 
(in our model) if the range is extended to include more of the lower 
main sequence (Table [4} . Indeed the change from the primordial 
value is less than 0.1. Baumgardt & Makino l l2003l) used TV-body 
models (not specifically geared to 47 Tuc) to assess the dependence 
of the mass-function index on mass lost by the cluster; and our re- 
sult is entirely consistent with their Fig. 9, since in our model 47 
Tuc has lost less than half of its initial mass. 

4.4 Pulsar accelerations 

FigdBshows that the model is very nearly consistent with observa- 
tions of the spin-down of pulsars in 47 Tuc. The fact that three lie 
just above the upper curve is consistent with estimates of the small 
intrinsic component (Sec |2.2.5l . The critical object in this plot is 
the pulsar wi th the largest nega tive period derivative. It is known 
as 47 Tuc S l lFreire et alll2003l) , and these authors showed that it 
implies a projected mass/light ratio M/L > 1.4 in the region of 
the pulsar. In fact the projected value of M/ L at the location of this 
pulsar in our model is about 1.1. Note, however, that our model is 
a little bright in the core (Sec |4.1b . which depresses the value. (In- 
cidentally, the projected value of M/L increases from this central 
value to about 2.3 at large radii; the global value is about 1.52.) 
Generally speaking, the tension between the requirements of the 
central surface brightness, on the one hand, and the acceleration of 
47 Tuc S, on the other, was the single most important constraint in 
our attempts to find a satisfactory model. 

4.5 Variation of Parameters 

Even though we had to work very hard to reach Model A, and dis- 
carded numerous models in the process, we can make no claim 
to the uniqueness of the model presented above, even within the 
limitations of the parameters which characterise our initial condi- 
tions. But in this subsection we wish to summarise the most notable 
changes (in the surface brightness, surface density, velocity disper- 
sion profile, mass functions and pulsar accelerations) which occur 
if we adjust our choice slightly. The list of parameters begins with 
those in Table [4] but we have also considered a number of other 
possibilities, which are included in the list below. Where we have 
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Figure 11. Comp arison between spindown rates of millisecond pulsars 
fcreire et al J2003h and the extreme contributions from the gravitational ac- 
celeration of Model A. 

varied a given parameter, the values of the remaining parameters 
are fixed and close to, but not always exactly the same as, those 
in Model A. Where no mention is made of changes, it means that 
the changes (in the five quantities listed above) are negligible in 
comparison with other entries. We generally desist from giving a 
theoretical interpretation of these empirical statements, each one of 
which might be a research project on its own. 

(i) N (Initial number of stars): an increase by 5% leads to an 
increase in the velocity dispersion by about 3%, and an increase in 
the maximum gravitational pulsar spinup by about 4%. 

(ii) n (Initial tidal radius): an increase by about 25% leads to a 
decrease in the central surface brightness by about 0.3 mag and a 
corresponding decrease in the luminosity function in the inner field, 
an increase in the outermost radius in the surface brightness and 
density profiles by about 10%, a decrease in the central line-of-sight 
velocity dispersion by about 5%, and a decrease in the maximum 
gravitational pulsar spinup by about 25%. Note that, because we 
hold the initial ratio of tidal and half-mass ratio constant, a variation 
of the initial value of rt also affects that of ru- 

(iii) rt /rn (Initial ratio of tidal and half-mass radii): An increase 
by 12.5% (i.e. a correspondingly more compact initial configura- 
tion, relative to the tidal radius) gives an increase in the central 
surface brightness by about 0.1 mag per square arcsec, an increase 
in the central velocity dispersion by about 4%, and an increase in 
the maximum gravitational pulsar spinup by about 18%. 

(iv) Wo (Initial King concentration): an increase by 1 leads to an 
increase in the maximum gravitational pulsar spinup by about 14%, 
but no other significant changes. This is also discussed in Sec |3.3l 
see especially Fig[3] 

( v ) fb (Initial binary fraction): an increase by 50% (e.g. from 
0.02 to 0.03) leads to a decrease in the central surface brightness 
by about 0.3 mag, a decrease in the central surface density by about 
19%, and a decrease in the maximum gravitational pulsar spinup by 
about 11%. 

(vi) «i (Index of the lower initial mass function): an increase 
by 0.1 (e.g. from 0.4 to 0.5) decreases the velocity dispersion pro- 
file by about 2%, and causes small changes in the outer luminosity 
function as expected from the discussion in Sec |3.3l (see especially 
FigH. 

(vii) Q2 (Index of the upper IMF): an increase by 0.2 (e.g. from 



2.8 to 3.0) leads to an increase in the central surface brightness 
by about 0.4 mag, an increase in the central surface density by 
about 0.15 dex, and an increase in the maximum gravitational pul- 
sar spinup by about 18%. 

(viii) nib (Break mass): an increase from 0.75 to O.85A/0 leads 
to a decrease in the central surface brightness by about 0.6 mag, a 
decrease in the central surface density by about 30%, some changes 
(in line with naive expectations from the change in break mass and 
the smaller central surface density) in the inner luminosity function, 
and a decrease in the maximum gravitational pulsar spinup by about 
20%. 

(ix) mi (Minimum initial mass): we did not experiment widely 
with this parameter, but can report that an increase from 0.08 to 
0.1 Mq causes an increase in the core velocity dispersion by about 
2%. 

(x) m2 (Maximum initial mass): a change from 50 to 150Mq 
leads to negligible changes. 

(xi) t (Time [Age]): an increase by lGyr leads to a decrease in 
the central surface brightness by about 0.15 mag, a decrease in the 
central surface density by about 20%, and a decrease in the velocity 
dispersion by about 6%. 

(xii) D (Distance): an increase by 0.5kpc leads to a correspond- 
ing reduction in the apparent length scale of the surface brightness, 
surface density and velocity dispersion profiles, and an increase in 
the inner luminosity function by about 8%. 

(xiii) Z (Abundance): an increase by 0.0005 from 0.003 to 
0.0035 leads only to a decrease in the maximum gravitational pul- 
sar spinup by about 4%. 

(xiv) <Tk (One-dimensional dispersion of natal kicks of neutron 
stars): a decrease from 190km/s (our canonical value) to 160 km/s 
leads to no significant changes, but a larger decrease to 130 km/s 
leads to an increase in the maximum gravitational pulsar spinup by 
about 10%. 



5 DISCUSSION OF THE MODEL 
5.1 The size of 47 Tuc 

A description of the spatial structure of globular clusters is often 
reduced to just three values: the core, half-mass and tidal radii. Our 
va lue for the tida l radius (Table |4) is rather comparable with those 
of lMevlanl ([1988, about 70 pc) and Table [2] (which converts to 56pc 
at the stated distance). Our value for the observational core radius 
(Table |4j is quite similar to that derived from Table [2] i.e. 0.52pc. 
There is, however, an amazing variation in quoted values for the 
half-mass radius. Data in Table[2]leads to a value of about 3.65 pc, 
and the difference from our value (Table|4j may be explicable if the 
former value should really be understood as a projected half-light 
radius. The values in Meylan's best models, however, lie close to 
9.5pc, and we have no explanation for such a wide disagreement. 
It does, however, have the effect of placing more mass at larger 
radii in his model, and this may be one reason why his fit with 
the velocity dispersion profile at large radii is more successful than 
ours (see Sec 14.21 ) 



5.2 Velocity anisotropy 

The velocity anisotropy of 47 Tuc near the cent r e has been 
constrained obser v ationa lly by iKing & Anderson! l l200lh and 
iMcLaughlin et al.l d2006h . The former authors also measured 
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Figure 12. Anisotropy in the proper motions of Model A and the observa- 
tional result for stars with my < 18.5 within 100" (McLaughlin et al. 
l200d) . The label on the ordinate specifies what is computed, the mean 
square radial and tangential proper motions at a particular line of sight 
being cr|j and Oq . 
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Figure 13. Projected fraction of main sequence binaries (as defined in the 
text) in Model A. The high points are an artefact of the way in which the 
surface density is computed, which can be very large if the line-of-sight 
nearly coincides with the radius of a binary. 



anisotropy at about 4' from the cluster centre, but gave no quanti- 
tative results. Fig|12lshows anisotropy among bright stars (V < 6) 
in our model, compared with the result for the brightest stars 
(my < 18.5) within about 100" of the centre ( IMcLaughlin et al.l 
l2006h . Though the observational result differs only marginally from 
zero, the results are consistent. A rapid rise at large radii, as seen in 
our mod el, is qualitative ly consistent with what was inferred indi- 
rectly bv lMevlar] dl988l) on the basis of model fitting. 



5.3 The binary frequency 

FigEl shows the radial distribution of the fraction of main se- 
quence binaries. This was calculated as in Sec |3.2.2l except that 
in calculating the surface density of single stars we include only 
main sequence stars, and in the surface density of binary stars we 
include only binaries in which both components are main sequence 
stars. In both cases also the apparent V magnitude was limited to 
the range 20.5 to 23.5. These are quite appropriate choices for com- 
parison with the observational data su mmarised in Sec |2.2.4| except 
that the range of magnitudes used bv lMilone et al.l d2008l) was a 3- 
magnitude range specified in J. In the radial range of their data, our 
result is just consistent with theirs. 

Of greater relevance for dynamics is the fraction of all bina- 
ries, not just those with main-sequence components. Our initial 
conditions generate a large proportion of relatively soft pairs. These 
are quickly destroyed, and the global fraction decreases by over 
60% (Table 0. Within the innermost 1000 stars, however, the bi- 
nary fraction increases with time as recorded in Figl 141 because of 
mass segregation. 



5.4 The dynamical role of dark remnants 

In the core, the velocity dispersion in the model is certainly af- 
fected by degenerate components, which make up about 5 1 % of the 
mass there, but their mass fraction decreases at larger radii to give 
a global ratio of about 34% (see Table[4 i. At face value, our resul t 
for the core contrasts with the finding of J 
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McLaughlin et alj l l2006l) . 



who suggested, on the basis of the profile of dispersion of proper 
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Figure 14. Fraction of all binaries in the innermost 1000 objects in Model 
A. 



motions, that the contribution of heavy remnants (such as neutron 
stars and white dwarfs) is only a fraction of a percent. The reso- 
lution of this discrepancy may simply be to understand that their 
estimate applied only to neutron stars, and perhaps the most mas- 
sive white dwarfs. Our global result is comparable to a number of 
ea rlier estimates, based on model fitting, which were summarised 
by iHeggie&Hutld 19961) . 

Despite our tinkering with the initial mass function, and 
despite giving natal kicks with a one-dimensional dispersion of 
190km/s, our model appears to produce a reasonable number (213) 
of neutron stars, . A recent indirect estimate, based on a quite 
simple demographic argument, p uts the total neutron star popu- 
lation in the cluster at "~300' ' dHeinke et al] d2005h: see also 
IVerbunt & Mevlanl ( fl988l) ). It has been argued dlvanova et al .l2008h 
that substantial numbers of neutron stars must be formed in pro- 
cesses leading to smaller natal kicks than assumed in our model, 
and so the number in the model could well be larger if this refine- 
ment were added. 

In our model there are only 19 stellar-mass black holes at 
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Figure 15. Evolution of lagrangian, tidal and core radii in Model A. Defi- 
nitions of the core radii are given in the text. 

the present day, and their total mass is small compared with that 
of other remnants. The maximum number present at early times 
is over 1000, but the number falls to less than 40 within the first 
30Myr. 

There is no sign of any earl y phase of runaway coalles- 
cence in the model. The results of IPortegies Zwart & McMillanl 
their eq.(17)), based on theory calibrated by TV-body mod- 
els up to TV = 64fc, imply that the collision rate is less than 
one per Myr. They also noted, however, that their result is likely 
to be an overestimate at half-mass relaxation times above about 
30Myr. Furthermore our upper mass function is steeper than theirs, 
and this helps to diminish the rate of collisions. Incidentally, 
the initial half-mass density, which is 0.28 x 10°, is within the 
range of a number of young massive clust ers at the present day 
jPortegies Zwart. McMillan. & GielesfcOld) . 

5.5 The dynamical status of 47 Tuc 

Because of its high concentration, it is tempting to assume that 47 
Tuc is close to core collapse. Though it is a massive cluster, its 
core relaxation time is given jHarris|[l996l) as t rc ~ 9 x 10 7 yr. 
For an equal-mass King model having the same concentration as 
47 Tuc at the present day, i. e. about 2.03, the core collapse time 
is about 400t re (from data i n lOuinlanlll996h . But it is known (e.g. 
IChernoff & Weinberg! [l990h that the time to core collapse in sys- 
tems with unequal mass is much smaller than that in equal-mass 
systems; these authors quote an example of a Wo = 7 King model, 
where a multi-mass system collapses in a time smaller than an 
equal-mass system by a factor of approximately 45. Putting all 
these factors together suggests that, indeed, the core of 47 Tuc 
should collapse within about lGyr. The large numerical factors in- 
volved, however, suggest that this argument is quite precarious. 

The spatial evolution of Model A is shown in Fig[T5] The core 
radius is represented by the two curves, one of them quite wavy, 
near the bottom of the figure, not far from the 1% Lagrangian ra- 
dius. These two curves are Bezier fits to rather noisy data on the 
core radius (cf. Fig|16l> generated in two ways. The curve labelled 
in the key as "dynamical" is determined from the mass-weighted 
three dimensional velocity dispersion and density of the innermost 
20 particles, using the formula r\ = 3{v 2 ) / (AivGpo)- The "obser- 
vational" curve is determined from a computation of the surface 



Figure 16. Evolution of half-mass and observational core radii in Model A 
and in a model with no primordial binaries but which is otherwise identical. 
The core radii are smoothed, but the original model data are also shown for 
Model A. Both runs were stopped at the point where the data ends. 



brightness profile, and is the radius at which this falls to one half of 
its central value (Sec l3.2.fl . 

We first notice that, except at the very start, there is a fairly 
steady (though decelerating) increase in the observational core ra- 
dius, at least in absolute terms. It also increases relative to the dy- 
namical core radius within the first few Gyr. We interpret the in- 
crease relative to the dynamical core radius as due to a combination 
of stellar evolution and mass segregation: as the turnoff mass de- 
creases, the observable core is dete rmined by star s of lower mass, 
which are more spatially dispersed. [Hurlev (2007) has described a 
similar effect in A^-body simulations, but attributed it to the core 
collapse of the heavy (and presumably dark) stars. Perhaps both 
effects are present, or perhaps they are simply somewhat different 
descriptions of the same processes. 

After the first few hundred Myr, the dynamical core radius 
is almost constant in absolute terms, but does appear to decrease 
slowly with respect to the half-mass radius (Fig|15t. The time scale 
of this process, which we refer to as core collapse, is difficult to 
estimate from Fig[T5] but the predicted future evolution is given in 
Fig[H] Evidently the collapse of the core slowly accelerates, but 
continues for at least a further 25Gyr! On this evidence, then, 47 
Tuc is far from core collapse. 

Now we turn to the mechanism which causes the overall ex- 
pansion of the system, the fact of which is evident from Figl 1 5 1 
Three processes come to mind: 

(i) mass loss from stellar evolution; 

(ii) primordial binaries; and 

(iii) new binaries, formed in three-body interactions. 

Empirical evidence comes from a comparison simulation in which 
the initial binary fraction was (Fig|16t; we refer to this as Model 
AO. The half-mass radius increases at almost the same rate in 
both models, even though Model AO contains at most one binary 
(formed in three-body interactions) until it reaches an age of about 
23Gyr (a few Gyr after the estimated time of core bounce, at about 
20Gyr). There are no binaries in it between about 500Myr and 
4Gyr, and yet the expansion of the half-mass radius of Model AO 
is very close to that of Model A, even during this period. It must 
be concluded, therefore, that the effect of binaries is almost negli- 
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gible in both models, and almost the entire expansion must be due 
to mass loss from stellar evolution, at least up to an age of 12Gyr. 

The difference between the two models is more noticeable in 
the core, however, and we conclude that the core of 47 Tuc would 
be about 20% smaller were it not for the action of the primordial 
binary population. In Model AO it seems clear that the core radius 
has passed its maximum by an age of 8Gyr, and the figure shows 
that core collapse is arrested at about 20Gyr. From then on three- 
body binaries form in greater abundance, and there is a slight but 
noticeable acceleration in the expansion of the half-mass radius. 

It may seem surprising that mass loss from stellar evolution is 
dominant. Its time scale, however, is set by astrophysics. The other 
possible competitors (core collapse, and heating by primordial bi- 
naries) scale with the relaxation time, other things being equal. This 
time scale increases with the mass of the cluster, and so for a suf- 
ficiently massive cluster, these stellar dynamical processes are rel- 
atively weak. 47 Tuc ra nks fifth amo ng the globular clusters of the 
Galaxy for luminosity dHarris|[T996t ). and, given the rather narrow 
range of M / L ratios, is certainly one of the most massive. 

It is an interesting and puzzling fact that the most lumi- 
nous Galactic globula r clusters tend to be the most concentrated 
dVan Den Bergr]|2003h . It is puzzling because, if we ignore stellar 
evolution, we expect the clusters of lowest mass to evolve fastest 
and therefore to approach or enter core collapse sooner; and this is 
the opposite of what is found. We see from this model of 47 Tuc, 
however, that high concentration is not necessarily any indication 
that the cluster is close to core collapse. 

5.6 Missing ingredients 

Several realistic aspects of star cluster evolution are omitted from 
our modelling. These include rotation, oblateness effects, a time- 
dependent tide (due to the ellipticity of the Galactic orbit of the 
cluster), tidal accelerations (as opposed to a cutoff), various aspects 
of binary dynamics (because we use cross sections for the inter- 
actions), and the complexities of initial conditions. Some of these 
omissions are common to virtually all simulations of star cluster 
dynamics, and have been since the field began. At the present day, 
however, especially in the context of globular star clusters, it is the 
question of initial conditions which is most controversial and most 
fluid. 

Initial conditions are complicated by a number of fac- 
tors, of which we highlight two. The first is early gas ex- 
pulsion. This has been t he subject of considerable research 
for some time now (e.g. Marks. Kroupa. & B aumgardt 20081; 
Baumgardt. Kroupa. & Parmentierl 2008 ; iGoodwin & Bastianl 
2006 ; iBastian & Goodwin! 120061 ; Fellhauer &Kroupal 120051 ; 
Gcvcr & Burked 1200 ll) . but we have ignored its effects. The 
second problem is that it has to be recognised that 47 Tuc cannot 
really be regarded as consisting of a single stellar population 
any more dAnderson et al.l l2009h - The appropriate paradigm for 
this situation, let alone its simulation, are full of uncertain ties 
dD'Ercole et ai]|2008l ; IChoi & Yil2008l ; iDowning & Sillsl2007l) . 

Faced with these complexities one might conclude that a 
modelling exercise like that in the present paper is meaning- 
less, misguided or at least premature. Our view on this is based 
on the fact that these missing features may affect only the first 
hundred million years of evolution, and perhaps less. We make 
the assumption that our models, after such a period of evolu- 
tion, do resemble hypothetical models which start with differ- 
ent i nitial conditions but do include the se missing ingredients. 
Thus iMarks. Kroupa. & Baumgardl d2008l) have shown that mod- 



Table 5. Initial central conditions for models of three clusters 



Cluster 



Central density Core radius Escape speed 





(M Q /pc 3 ) 


(pc) 


(km/s) 


M4 
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0.31 
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3.0 x 10 6 


0.22 
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47 Tuc 


1.0 x 10 6 


0.37 
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els which include gas expulsion produce a mass function which ap- 
proximately resembles what we have chosen for our IMF (at least 
qualitatively, in the sense that the lower mass function will have an 
index smaller than the canonical value of around 1.3). 

These processes are likely to depend on the mass of the cluster, 
and it is not surprising that our choice of IMF for 47 Tuc is quan- 
titatively somewhat diff erent from those we preferred for the clus- 
ters N GC6397 and M4 dHeggie & Gierszll2008al ; lGiersz & Heggie] 
120091) , which we consider to have been less massive initially by a 
factor of about four. But the length scale is also important. It is in- 
teresting to note, for instance, that the central density of our initial 
conditions for these three clusters lies within a fairly small range 
(Table[5}. The same is true of the central escape speed. Indeed, for 
most of its life so far, the central escape speed in 47 Tuc is larger 
than for the other two clusters by a factor of only about two (see 
iGiersz & Heggiell20od Fig. 12), though it still exceeds 40km/s at 
the present day. If chemically contaminated winds are a significant 
factor in the production of a second generation, the retention factors 
among different clusters may not vary over a wide range. 

When this project of fitting Monte Carlo models to individual 
globular clusters was begun, it was thought that post-collapse clus- 
ters would be most difficult to treat, because of the complex inter- 
action of several processes which would lead to their present struc- 
ture. Our conclusion now is that pre-collapse clusters are equally 
difficult, because they have changed less from their virtually un- 
known initial conditions. 



6 CONCLUSIONS 

We have constructed a dynamical evolutionary model of the mas- 
sive Galactic globular cluster 47 Tucanae (NGC 104). The model 
takes into account dynamical interactions between stars and bina- 
ries, the stellar evolution of these components, and the effect of the 
Galactic tide. We make no claim for the uniqueness of this model, 
though numerous other models were eliminated in the search for it. 

The model begins with a moderately concentrated King model 
(Wo = 7.5) without primordial mass segregation. The initial and 
present-day characterisation of the model, including values for its 
mass, spatial dimensions, the mass of degenerate remnants, and the 
binary fraction, are given in Table [4] In summary, the initial mass, 
half-mass radius and number of stars are 1.64 x 10 6 M Q , 1.91pc 
and 2.00 x 10 6 , respectively. By the present day these values have 
changed to 0.90 x 10 6 M Q ,4.96pc and 1.85 x 10 6 , respectively. 
The most unorthodox aspect of the initial conditions is the mass 
function, which is steeper than Salpeter for high masses, and rela- 
tively flat at low masses. 

We judge the success of the model by comparison with obser- 
vations of 47 Tuc at the present day: the surface brightness and sur- 
face density profiles, the velocity dispersion profile, the luminosity 
function at two radii, and observed pulsar accelerations. The centre 
of the model is a little bright, perhaps because of some flaw in the 
recipes for stellar evolution of post-main sequence stars. The least 
satisfactory comparison is with the velocity dispersion profile. On 
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the other hand it is not clear how to assess the velocity dispersion 
profile of a spherical, non-rotating model which is constructed so as 
to resemble the spatial structure of a somewhat flattened, rotating 
object. 

We find that the primordial binary population of 47 Tuc is 
playing a significant role in the evolution of the core radius but not, 
at present, a dominant one. The cluster is far from core collapse, 
which will not take place for another 25 Gyr at least. 
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